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SUMMARY 


A method for computing one -dimensional, steady, nonequilihri'um channel 
flov of an arbitrary mixture of diatomic and monatomic gases is presented. 

The method assumes the chemical and vibrational nonequilibrium rate processes 
are not coupled. The nimerical formulation permits accurate evaluation of all 
variables appearing in the rate equations, and it precludes the necessity for 
computing expressions vhich become indeterminate in the near -equilibrium 
region. 


INTRODUCTION 


The flow of a nuilticomponent gaseous mixture which is influenced by 
thermochemical nonequilibrium has proven to be difficult to predict in the 
near -equilibrium region. In the analysis of a flow of this kind, the equa- 
tions which describe the fluid dynamical processes must be solved simulta- 
neously with equations which describe the nonequilibrium rate processes. 

These rate equations take the form of a coupled system of nonlinear differen- 
tial equations, and in order to proceed toward a prediction of the flow, it is 
necessary to use n'umerical methods. 

A numerical solution of the equations that govern the flow is possible in 
principle, but has not been demonstrated to be practical for all cases that 
may arise. Early investigators found excessive amounts of computation neces- 
sary in order to achieve numerical stability in the near -equilibrium regions 
(refs. 1, 2 , 3 ). The instability arises in part from the necessity for eval- 
uating expressions which, as equilibrium is approached, become indeterminate 
in the sense that their straightforward calculation leads to large numerical 
errors. Several methods (refs. 4, 5^ 6 ) that have been of practical value 
were devised for the approximate numerical integration of such so-called 
"stiff equations, " although they did not entirely remove the underlying cause 
of the difficulties. In these early studies, the interest was centered pri- 
marily on the physics of the problem, and the n'umerical difficulties were 
regarded as an inconvenience. 

Recently, three new methods have appeared that greatly reduce the comput- 
ing time in the near-equilibrium regions of the flow (refs. 7^ 9)* These 

later papers indicate a growing interest in the mathematical problem which, 
evidently, has not been solved satisfactorily. 



The nature of the problem may be defined by recalling that regions requir- 
ing special treatment are frequently encountered in numerical calculations. 
These are usually regions in which some physical quantity approaches an asymp- 
totic value. A standard remedy is to find a closed -form solution of a simpli- 
fied set of equations that apply in the asymptotic region. A difficulty at 
first appears to preclude the use of this device in the present problem. The 
trouble arises from the fact that the equations are nonlinear over extensive 
near-equilibriirm regions so that a closed-form solution cannot be found. The 
method, however, can be applied if the equations eire linearized independently 
in each small step of the numerical calculation, as was done by Moretti 
(ref. 8). At first, it appears unwieldy to Incorporate such a procedure in a 
numerical integration program, but the resulting reduction in computing time 
gives good reason to do so (refs. 8, 9, 10). Also, such a linearization pro- 
cedure provides the only method available at present for computing accurately 
the so-called reaction variables that appear in the rate equations. 

The present method differs from that of Moretti (ref. 8) in that accurate 
evaluation of the reaction variables is included. To accomplish this, consider- 
able manipulation of the original rate equations is necessary. However, upon 
completion of these steps, a system of equations is obtained that applies for 
equilibrium, near -equilibrium, or nonequilibrium flows. The method to be pre- 
sented is based on the requirement that means be found for computing, as accu- 
rately as might be specified, all quantities appearing in the equations. This 
requirement was previously employed in conjunction with a simplified model for 
air (ref. 10) and has been discussed in connection with more general reacting 
systems (ref. 11). 

Specifically, this report contains a method for computing one -dimensional, 
steady, nonequilibrium channel flow. The method pertains to an arbitrary 
multicomponent mixture of ideal diatomic and monatomic gases, free from trans- 
port effects and from interaction between the chemical and vibrational non- 
equilibrium rate processes. The numerical formulation mahes unnecessary the 
computation of expressions which become numerically indeterminate in the near- 
equilibrium region, but it requires the solution of an eigenvalue problem at 
each field point of the one -dimensional flow. 

The numerical difficulties inherent to calculating the eigenvalues of sin 
arbitrary matrix are well known. However, with the numerical formulation of 
this report, the matrix in question is simileir, at every field point, to a 
real symmetric matrix which is easily constructed. This fact is used to 
greatly simplify the computational problem. 


SmBOLS 


A matrix of coefficients (see eq. (A1)) 

A“^ inverse of matrix A 

A^ transpose of matrix A 
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A(x) cross-sectional area of the channel at point x 

D(t) diagonal matrix with diagonal elements e (i = 1 ^ 2 ^ . . r) 

diagonal matrix with diagonal elements ^ (i = 1^2, • . .^r) 

E matrix whose column vectors are the eigenvectors of matrix A 

energy of dissociation per mole of species J 
internal energy per mole of species j 

ij combined rotational^ electronic, and translational energies per mole 

of species j 

e^ vibrational energy per mole of a given species 

e^{T) equilibrium vibrational energy per mole of a given species at 
temperature T 

e^j vibrational energy per mole of species j 

f natural vibrational frequency for a given molecular species assiiming 

it is a classical harmonic oscillator 

H matrix (see eq. (38c)) 

h enthalpy per unit mass of the gas mixture 

hxj element in the ith row and jth column of the matrix H 

(see eq. (38c)) 

ho Planck* s constant 

I identity matrix 

k Boltzmann* s constant 

^bi backward rate coefficient of the ith reaction 

kfi forward rate coefficient of the ith reaction 

Kj_(T) equilibrium constant of the ith reaction 

kn n+i probability that a molecule vibrating in the nth energy level will 
^ make a transition to the (n + l)th level in unit time 

M matrix (see eq. (38b)) 


chemical formula of species j 

element in the ith row and jth column of matrix M 
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Nq Avogadro * s nimber 

^kj total niiniber of atoms of element k in species j 

p pressure 

Qij rate of production of species Mj by the ith reaction per unit 

volume 

R universal gas constant 

r number of chemical reactions 

s number of chemical species 

T temperature 

t time 

U see equation (27) and explanation thereafter 

u speed of flow 

V see equation (27) and explanation thereafter 

X distance along the axis of the channel 

aj_ constant appearing in equation {kO) 

Pi integer equal to vi^ - Vi 

Pij integer equal to v±^ - Vij 

7j concentration of the Jth species in moles per unit mass of the gas 

mixture 

Aeyj increment in eyj during the interval At 

AT increment in T during the interval At 

At incremental interval of time 

A7. increment in y a during the interval At 

J ^ 

Ap increment in p during the interval At 

AXi increment in Xj_ during the interval At 

€n energy of the nth vibrational level of a molecule 

A diagonal matrix whose elements are the eigenvalues of the matrix A 


4 



M-n 

Vi 


Vi* 


ith eigenvalue of the matrix A 

number of molecules of a given species which are in the nth vibra- 
tional energy level 

s 

integer equal to 




s 


integer equal to 








v! . 


P 

0 




^i 

-► 

X 


stoichiometric coefficient of species Mj on the left side of chem- 
ical equation i 

stoichiometric coefficient of species on the right side of chem- 

ical equation i 

density of gas mixture 

sum of all species concentration coefficients^ 7 j (j = 1^ 2^ . . s) 
vibrational relaxation time for the jth species 
reaction variable (see eq. (5)) 

column vector whose elements are Xj_ (i = 1^ 2^ . • r) 


Subscript 

( )q used to indicate the numerical value of the quantity taken at the 
beginning of the interval At 


ANALYSIS 


Although they have been developed elsewhere, the basic equations for non- 
equilibrium channel flow will be derived here. This derivation will provide 
a background for the numerical method used in their solution. 


Chemical Rate Equations 


It is assumed that 
involved in r chemical 


s chemical species, (j = 2, . . ., s) are 

equations. These equations are represented by 
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kbi 


V • .M . 

ij J 


( i — 1^ 2f . . . ^ r) 


J=i 


( 1 ) 


where Vij and are stoichiometric coefficients, the Mj are the chemical 

formulas of the species, and kfi and kbi are the forward and backward rate 
coefficients, respectively, of the ith reaction. 


The ith reaction produces a quantity of species Mj at the rate Qij 

(moles/unit volume/unit time). For a fluid without gradients of concentration, 
pressure, or temperature, the law of mass action leads to the following 
expression for Qij 





ij 




( 2 ) 


In this equation, 7j (moles/unit mass) is the molar concentration of species 
Mj ( j = 1, 2, . . . , s) , p is the fluid density, and 


V 


i 



s 



0=1 


and = vlj - Vij 


It is asstimed that equation ( 2 ) holds for a flowing gas. 

The equilibriijm constant for concentrations of the ith reaction, Ki(T) 
is a known function of the temperature, T, and at equilibrium. 


Ki(T) 


kfi(T) 

kbi(T) 


(i = 1. 2, 


• ^ 


r) 


It will be assTomed that this relationship holds when the reaction is not in 
equilibritim as well. 

The total rate of change of 7^ is given by 

J 

d 7 j 

p ^ =" 2^ %j (j = 1, 2, . . ., s) 
i=i 
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Thus 


d7J 

P aF 


I ^ 4 ^/ 


^ ^11 


7i 

7=1 


i=l 


V 




! s V 
^ Tl 7i 

7=1 ^ 



(J = 1, 2, 


s) 


(3) 


These equations can he written 



1=1 


(j = 1 , 2 , 


s) ik) 


vhere the dimensionless variables Xj_ are defined by 


Xi - 1 


1 

Ki(T) 


r S 317 

7z 

Z=^i 


(1 = 1 , 2 , 


r) 


(5) 


and 3i = - Vj^. (See ref. 3*) In near-equilibrium regions, the 

approach zero, aind they cannot be computed accurately from equations ( 5 )* 

Their accurate evaluation is essential if equations (4) are to be used, since 
the d 7 j/dt cannot be set equal to zero in near-equilibrium regions. In the 

method to be developed here, equations (4) will not be evaluated numerically, 
but the method is capable of producing accurate values of the Xi if they are 

desired. When equations ( 5 ) are differentiated with respect to time, there is 

obtained a set of relations between the Xj^ (i = 1 , 2 , . . ., r), the 7 . 

(j = 1, 2, . . ., s), and the fluid dynataics variables p and T. 


dXi 

dt 


1 Pi 

kT^tT p 


S 



Z -1 




1 

Ki^(T) 


(- 


dT J 




s Pi7\ 
n 7i ) 
7=1 / 


1 

Ki(T) 


s Pl7/Pi\ 
■ n 7l (— ) 
7=1 ^ / 



(i = 1 , 2 , . . r) 


(6) 


Equations ( 6 ) are useful because they contain no forms which become indeter- 
minate at some point in the flow field, including the near -equilibrium region. 


T 



Equations (4), (5), and (6) can be confined to obtain another set of 
equations which will also be used in the computation of the flow. 


V ( 

^ ^ aijXz + F(T,Xi) ^ + G(p,Xi) ^ (i = 1, 2, . . . , r) (?) 


where 






J=i 


and 


1 = 1, 2, . 

2 = 1, 2, . 


r 

r 


F(T,Xi) = (1 - Xi) 


Ki(t)J V dT 


dK^ 


7 ^) (i = 1, 2, . . r) 


<^(p^Xi) = (^i “ l)f “d” 


(i = 1, 2, , . r) 


Vibrational Rate Equations 

It is assumed that the molecules of the gaseous mixture are free to 
vibrate only at a discrete set of energy levels 

€n = (^n f l^hof (n = 0, 1, 2, . . .) 

where h^ is Planck* s constant and f is the natural frequency for the 
species in question assmning it is a classical harmonic oscillator. It is 
further ass'umed that changes in energy can occur only between adjacent levels. 

The probability that a molecule vibrating in the nth level will make a 
transition to the (n + l)th level in a tinit time is designated by k^ n+i‘ 
With this notation, ko^i and kj_ q are the probabilities of excitation^ and 
de-excitation per unit time, respectively, between the ground state and the 
first excited energy level. Let ij.^ the number of molecules in the nth 
vibrational energy level. Then 


<i^^n 

dt 


kn-i,nlJ-n-i ” ^n,n-iM’n + kn+i,nM’n+i 


^n,n+i^*’n 


( 8 ) 
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In reference 12 it is developed that 


1 


= exp 



( 9 ) 


where k is the Boltzmann constant. For other energy levels 


kn,n+i 


= ^n+i,n 



and from the quantum-mechanical study of transition probabilities 


^n^n-i - 


( 10 ) 


If equation (8) is combined with equations (9) and (lO), it is possible to 
oTotain 


diJ.. 

dt 


In r 

~ = + l)l^n+i " ^^^n + 


exp 


-hof 

kT 


I [n|in-i - (n + l)nn]j" 


(11) 


Let Nq be Avogadro’s number^ and 


00 

w = ^ Hn 

n=o 


If is used to denote the vibrational energy above the ground (n - o) 

state of a given species^ per mole of species 


^ V 

^v ]\j- ) ^nM'n 

n^o 


Nohof ^ 

N L 

n=o 




Therefore 


00 

de^ ^0^0^ ^ 

^ ^ dt 

n^i 

VThen equation (ll) is combined with this last equation^ it is possible to 
manipulate the resulting expression into the form 
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Nollof* 


The expression 




is defined to be the vibrational energy^ per mole of species^ corresponding to 
equilibrium conditions at temperature T. Also^ the relaxation time^ is 
defined by the relation 


- - k. 


1 - exp 


With these definitions^ equation (l2) can be written in a more compact form. 

de^r [^v(t) - 6v] 


Since k^^o is a function of density^ temperature^ and chemical species^ t 
is a function of these variables. It is customary to write^ for the jth 
chemical species^ the vibrational nonequilibrium ecquations in the form 


[e-yj(T) - Pyj] 
Tj(p,T) 


(j = 1, 2, . 


The vihrational energy of species J is related to the total internal energy 
per mole, ej, of species j by the equation 




vhere ej(T) is a function of temperature which expresses the rotational, 
electronic, and translational energies per mole of species j. The term 
a constant, is the energy of dissociation per mole of species j. 
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Fluid Dynamical Equations 

Since the flow under consideration is to he one -dimensional and steady, 
all quantities used in its description are functions of a single independent 
•variahle. The channel is presumed to he described hy some function of x, 
the distance along the axis of the channel from some fixed point. However, 
the flow variables will he considered as functions of time, t. In the conser- 
vation equations which follow, the speed of the flow is represented hy u, the 
pressure hy p, the cross-sectional area of the channel hy and the enthalpy 
per unit mass hy h. The conservation equations are 


Mass 

puA = Cl 

( 1 ?) 

Energy 

1 p 

- U*^ + h = C2 

(18) 

Moment-um 

1 2 r dp 

2 U + / - = C3 

(19) 


The numbers ci, c^, and C3 are constants of integration. 

It is assimed that either a function A(x) or a function p(x) is given. 
When A(x) is given, the cross-sectional area of the channel is prescribed; 
when p(x) is given, the pressure is a known function along the channel. 


For a mixture of gases composed of s species, 
state is 



the thermal equation of 


(20) 


In this equation, R is the -universal gas constant. In later equations, the 
notation ^ 

I ’■3 ' ° 

j-1 


will he used. Thus 


p = ap R T 


(21) 


By differentiating equation ( 2 l) with respect to time, and using 


^ d_ _ 
dt dt dx ~ ^ dx 
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it is found that 


s 


^ ^ T ^ 

P dx ~ dt P dt 


a 


0=1 


dt 


(22) 


Equation (22) is useful when p(x) is specified. When A(x) is given 
instead, equation (22) can he adapted hy using the equations for conservation 
of mass and moment-urn. To this end^ the conservation of mass equation can he 
written 


(uA) ^ + (pA) ^ + (pu) H = 0 (23) 

and from the conservation of momentum 


or 


u 


du ^ 1 dp 
dx p dx 


0 


^ = 1 ^ 
dt “ P dx 


(24) 


When equations (23) and (24) are combined, there is obtained 


dp dp P'u-^ dA 

dx “ ^ dt A dx 


(25) 


Finally, by substituting the right side of equation (25) for dp/dx in equa- 
tion (22), it is seen that 


s 

u3 ^ ^ M T / uf£-\ dp T y ^'^3 

ARC dx dt P V P y dt a ^ dt 

j=i 


(26) 


This equation will be useful when A(x) is specified. 
In general 


^ T ^ T 

dt P dt a 



= V 


(27) 


When A(x) is specified, U = [1 - (u^/oRT) ] and V =• (u®/ARc ) (dA/dx) . When 
p(x) is specified, U = 1 and V =• (uT/p) (dp/dx) . 
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Eq^uation (2j) is one of the two fluid dynamical equations which will be 
used in the numerical solution of the nonequilibrium flow. The other such 
relation is derived from the equations for conservation of energy and momentimi. 


The specific enthalpy, h, is given by 


h = e + I 

For a mixture of s chemical species, the internal energy, e, per unit mass 
is 


s 





Therefore 


s 





(28) 


When this equation is differentiated with respect to time, there is obtained 


^ V ^ V CRT 

dt ~ ’'j dt i dt P dt " P dt 

j=i 


(29) 


When the equations for the conservation of energy and momentum axe differen- 
tiated and then combined, it is possible to show that 


to ^ 1 ^ 
dt P dt 


(30) 


From equation (l6). 


_ <3.ej 

dt dt dt 


(j = 1, 2, 


s) 


and since dej/dT is a known function of T, it is convenient to write 


dt " dT dt dt 


(j = 1, 2, 


• y 


s) 


(31) 
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When equations ( 29 ); (30) ^ (31) axe combined, there results 



de 


VO 


dt 


a 

I 

0=1 


d7^ 

t 

'j dt 


oftT dp 
~ dt 


(32) 


Equation ( 32 ) together with equation ( 27 ) will be used with the equations of 
the nonequilibrium rate processes in the numerical procedure which is to be 
described. 


MJMERICAL PROCEDURE 


For one -dimensional channel flow, the field points can be stationed along 
a line which extends the length of the channel. The variables p, T, /j, and 
Svo (j = 1^ 2, . . ., s) from which all other flow variables can be obtained 
must be computed at each of these stations. 

As a fluid element moves downstream from one station to the next, an 
interval of time At elapses, and there axe corresponding increments Ap, AT, 
A 7 j, and Aeyj. Let it be assumed that the fluid element is at a given station 
at time t© and at the next downstream station at time (to + At). Then 


p(to + At) = p(to) + Ap 


T(to + At) = T(to) + AT 


7j(to + At) = 7j(to) + A7 j 
® vj('to + At) = e-vj(to) + Ae-vj 


(j — L, 2, . . ., s) 
( j = L, 2, . . • , s ) 


The computational problem is to find numerical values for the increments Ap, 
AT, A 7 j, and Aeyj when numerical values of the flow variables axe known at 
the upstream station. Since values for all quantities axe assiimed known at 
the beginning of the interval (to^ to + At), there is no loss in generality- 
in taking to = 0 and denoting values at the beginning of the interval by the 
subscript o. Thus 

p(At) = po + Ap 

T(At) =. To -t AT 


7j(At) = (7j)o +A7 j (j - 1, 2, 


s) 


®vj (At ) — ( ®vj ) o ^®vj ( j = L, 2, 


s) 
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Computation of the Increments Aeyj 

As a consequence of the assumption that vibrational and chemical relaxa- 
tion rates are not coupled^ an expression for the Aeyj . . .^s) 

can be found independently of the increments ^ 7 jo Since p and T are func- 
tions of time^ equations (15) can be -written 



(eyj - evj)p 


(j = 1 . 2 , 


s) 


It might be presumed^ then^ that the increment Aeyj 

(Svj - ®vj)o 


Ae' 




(T.i) 




At 


(j = 1, 2, 


could be computed from 
» • • } ® ) 


Hovever_, use of these equations are not practical because of the behavior of 
the right-hand side in the near -equilibrium region of the flow. In the near- 
equilibrium region^ very large while (eyj - eyj) is very small- This 

makes an accurate evaluation of Aoyj impossible even for small At. Apart 
from this^ the use of " ®vj)o over the entire interval^ At^ is of ques- 

tionable validity. The quantities Tj^ Oyj^ and Oyj can be described as 
slowly varying in the near -equilibrium region^ because their change during a 
short-time interval is small enough that their value during the interval can 
be approximated by a constant. On the other hand, (Cyj - eyj) can change by 
several times its own magnitude during a small interval of time when the 
upstream station is in the near -equilibrium region. 


A means of avoiding these numerical difficulties is provided by the 
behavior of Cy j . Since it is a slowly varying function, de’yj/dt can be con- 
sidered constant over an interval At. Likewise, Tj can be considered con- 
stant over the interval. 


With this in mind, the following equations are derived from equations (l 5 )* 


de. 


vj (ievj 


dt 


dt 




(j 2 , . . ., s) 


These can also be written 


it ^ ^ 


1 




i) + 


de 




dt 


2 , . . s) (33) 


Over the interval At, Tj = (”^0)0 


de 




Ae 




dt 


At 
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Hence^ oyer this interval^ equations (33) hecome linear differential equations 
with constant coefficients. Their solution is found to be 


Ae 


[evj(t) - evj(t)] = “ ®vj) 


At 




Aeyj 

At 




(j-1^ Zf • • s) 

where 0 ^ t ^ At. From this solution, the Increment in (eYj - ©vj ) corre- 
sponding to At can be found 





At J 


(j = 1, 2, 


,, s) 
(34a) 


Equations (S4a) show that as (tj)q “*■ 0 (equilibrium is approached), the 
increments A(eYj - eYj) -*■ - (©vj “ ®vj)Q* There is no difficulty in computing 
these increments in the near -equilibrium region by means of equations (s4a) . 

The increments Asyj can now be determined accurately. 


AeYj - AeYj 


and 


AeYj - 


- A(e 
‘^®vj\ /^AT 


dT 


At 


'vj ■ =vj) 
deYj' 


( j — 1, 2, . . . , s ) 


At = 


dT 


AT 


(j — 1, 2, . . ., s) 


These last two sets of equations are combined with equations (S4a) to give 



(j = 1, 2, . . ., s) (34b) 
Computation of the Increments 

When the increments AeYj (o = 2, . . ., s) have been found by the 

method just outlined, there remain the (s + 2) increments Ap, AT, and A/j 
(j = 1, 2, . . s) yet to be computed. These will be found by solving a 
system of (s + 2) linear equations in which the variables are the (s + 2) 
increments in question. However, these equations contain the increments AXq 
(i = 1, 2, . . ., r) taken by the variables Xq during the interval At. It 
will be shown that upon solution of an eigenvalue problem, each AXq becomes 
a linear function of Ap and AT. 
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The variables Xx 
equations (7) 


dXi 

at“ 


(i = 1, 2., . . r) satisfy the coupled system of rate 
r 

- ^ aijXz + F(T,Xi) ^ + G(p,^i) ^ 


where 

^7 = (l - Xi) 




-X B Vta 

n 7 .^^ 

j=x J 


/i = 1, 2, 

V = 1. 2, 



F(T,Xi) 


= (1 - Xi) 


1 

Ki(T) 



G(p,Xi) = -(l - Xi)l 


Pi 


¥ith the exception of itself^ all the variables which appear in the coef- 

ficients a^x^ F(T,Xi)^ and G-(p^Xj[) are assimed to he slowly varying in the 
sense that they change hy only a small percentage of their magnitude during 
an interval At. Hence^ they will he approximated hy a constant value over 
this interval. In the near-equilihrium region^ however^ the variables Xi 
(i = 1^ 2^ . . r) can change hy several times their own magnitude during 
the interval. This magnitude must he small because Xi 0 as equilibrium is 
approached. However^ in any region of channel flow^ including the near- 
equilibrium region, (l - Xi) can change hy only a small percentage of its value 
during a small interval of time. This fact allows the coupled^ first-order 
system, equations (7), to he "linearized” hy talcing the coefficients to he con- 
stant over the interval At. Equations ( 7 ) are then written 


dX 


- = - £ (aii)^Xz + [F(T,Xi)]o(^^) + [G(p,Xl)]o ^ 

(i = 1, 2, . . r) (0 < t <At) ( 35 ) 


Z-1 


When the constants Ni are defined through the equation 

[F(T,Xi)]^ [G(p,Xi)] 


Ni - 


AT + 


Ap 


(i =. 1, 2, . . r) 


At • At 

the differential equations (35) can he written in vector form 

dX 


dt 


= -AX + N 


(36) 
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where 



|X,1 




/(aii)o 

(^12)q • • • 

(air)o 


' X2 


[ 

N2 


(^21)0 

(^2)0 * • • 

(^2r)o 

X 

• 

, IT = 

• 

y and A = 

• 

• • • • 

• ■ • • 

• 







(^r2^Q • ' • 

(arr)o/ 


Although numerical values cannot he assigned to Ni (i - 1 ^ 2 ^ . . r) until 
AT and Ap have been founds this must not obscure the fact that each Nj_ can 
be considered to be a constant over the interval At. 

If x(t) is a vector flection which satisfies the differential equa- 
tion (36) such that Xo = then each (l - 2 ^ . . r) is deter- 

mined by the equation 


^ = X(At) - X( 0 ) 

Consequently^ a general solution must be found for equation (36). 

Toward this end^ a transformation will be used to replace equation (36) 
with an uncoupled system. Theorem 3 the appendix shows that such a trans- 
formation exists. This theorem states that there are eigenvectors^ r\^ 

(i = l^ 2 ^ . . .y r) y of the matrix A which are linearly independent at every 
point in the flow. A simplified method for finding the eigenvalues and eigen- 
vectors of A is provided by theorems 1 and 2 in the appendix. 

If E is the matrix whose ith column is the eigenvector which 

corresponds to the eigenvalue 7\± of Ay then E is nonsingular since its 
column vectors are linearly independent (see ref. 13 )* Moreover^ 


A = E A E“^ 


where A is the diagonal matrix whose diagonal elements are the eigenvalues 
of A. 


A - diag (Ai^ 


Ai 


• ^ Ar) = 


A2 


Ar^ 


18 



Equation (36) can be written 


dt 


-(E a E“^)X + N 


or 


E"^ = -AE"^ + E~^ 

dt 


Setting ^ = 


E X yields 


^ = -A? + E"^ 

dt 


( 37 ) 


Equation (37) is an uncoupled system because A is a diagonal matrix of 
constants. 

If \l;(t) is a general solution of equation (37) then 

X(t) = E\|/(t) 

is a general solution of equation (36). 

^ The general solution^ ^(t)^ is found by obtaining a general solution 
^h.(t) of the homogeneous equation 


dt 




then adding to it any particular solution, 1^p(t) of the equation (3T)» 




The solution, t^(t), of the homogeneous equation is easily verified to be 


where 


^l^(t) = D(t)C 


c = 

fS 

and D(t) - 

e-Aat 


\:vj 
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A particular solution, ?p(t) of equation (37) ^ can be obtained by solving 

0 = -Atp + E"^N 

This gives 

?p = Di(t)E‘^N 

where the diagonal matrix Di can take two possible forms. If all the 
eigenvalues of A are nonzero 

Di(t) = A-i = diag ^ ^ 

The possibility that some of the eigenvalues may vanish, however, cannot 
be excluded. This will happen, for example, when there are more chemical 
equations than distinct chemical species. Each eigenvalue of A which 
vanishes must have its reciprocal in Di(t) replaced by t. For example, if 
the second and rth eigenvalues vanish, Di(t) takes the form 

D,(t) .aiag(i,t, 

Finally, 

?(t) = ^p('t) 

= D(t)C + Di(t)E“^N 


Therefore 


x(t) = [ED(t)]C + [EDi(t)E‘^]N 

It is possible to determine the constants Ci^ from the 

equation 

Xq = x(0) = [ED(0)]C + [EDi(0)E"^]N 

From the definition of D(t), it is seen that D(0) = I, the identity matrix. 
Consecquently^ 

C = E‘% - [Di(0)E"i]N 

and this yields 

X(t) = ED(t){E"% - [Dx(0)E"^]f} + [EDi(t)E"^]N ( 38 a) 

The increments, (i = 1, 2, . . ., r), for the interval At are found 

from 


^ = X(At) - 
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Therefore 


^ = [ED(At)E ^ - I]Xq + [EDi(At)E~^ - ED(At)Di(0)E‘^]N 

For a fixed value of At^ the diagonal matrices D(At) and Di(At) become 
matrices vith constant elements. For purposes of notation^ let 


[ED(At)E~ - I] M = m^j (i^ j 1^ 2^ . . . ^ r) (38b) 


and 


[EDi(At)E“^ - ED(At)Di(0)E"^] = H = hij (i, 0 = 1, 2, . . r) 


( 38 c) 


Then 


AX = MXq + HW 


In terms of the increments Ap and AT, this is written 


^i = ^ + 

0=1 


h. 




■ij 


At 


Lj=i 


AT + 


h. 


P^X-i ) 




IJ At 


>-0=1 


Ap 


(i — 1, 2, . . r) ( 39 ) 


Computation of the Increments Ap, AT, and A/j 
With equations (6) and (39) combined there is obtained 


1 dKj 

Ki(t) dT 


OC-i 


y f(t,x^)< 




At 


AT - 


Pi 


+ a± 


G(p.^j)o 

hiZ 


-Z=ii 


Ap 


where 


] ^7j - ^i ^ ^iz(X7^)o (i - . • y 


J = i 


Z=.i 


ai - 


Kq(T) 


- J=i ^ 


(1*0) 
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The two linear equations below are easily obtained from equation (27) and a 
combination of equations ( 32 ) and 

s 

AT + A7j + (v Ap = (v)q At (la) 

° J=l 



s s 

+ Y (®j)o j^. “ 4 ~ 

J=1 j=l 

ih-2) 

Equations (40), (4l), and (42) comprise a system of (r + 2) linear equations 
in the (s + 2) unknowns AT, Ap, and A 7 j ( J = 1, 2, . . . , s) . If it were 
possible that r = s and that equations (40) form a linearly independent set, 
these (s + 2) equations can be solved for the increments in question. However, 
both conditions cannot be satisfied simultaneously^ and supplementary linear 
equations in the variables AT^ Ap^ and A/j ( j ^ 1 ^ 2 ^ . . s) are needed. 
These supplementary equations can be obtained^ just as in equilibrium from 
atom conservation equations. 

The system of equations (4-0) always contains a linearly independent subset 
which can be determined by reference to the original chemical equations ( l) . 

The chemical equations (l) can be ordered so that each succeeding equation^ 
after the firsts contains a chemical species not appearing in any of the 
previous equations. If this ordering exhausts all of equations (l)^ there will 
be some number^ say ri (where ri < r) ^ of equations (40) which is linearly 
independent. This linearly independent subset of equations can be determined 
before any computation is begun^ and the required atom conservation equations 
can be adjoined. 

In the s species^ there are a fixed number^ m^ of different elements. 

Let n]£j be the total number of atoms of element k (k = 1^ 2^ . . m) 
in species j. The total number of atoms of element k must be constant so 
that 

s 

Z ^7 j 

n-kj ^ 0 (k 1 , 2 , . . . , m) 

j=i 

When d 7 ./dt is replaced by A/j/At and both sides of this equation are 
multiplied by At^ there is obtained 
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s 




(k = 1, 


I 




2^ • • • ^ m) 


From these m linear equations^ the required supplementary equations can he 
obtained. 

It should he noted that equations (l+O) were derived independently from 
the fluid dynamical equations (l?)^ (l8), and (l9)* They can he used^ there- 
fore^ in conjunction with other sets of fluid equations which arise in, for 
example, three-dimensional flow, or unsteady flow. 


CONCLUDING REMAEKS 


In order to compute accurately the nonequilihrium channel flow of a 
multicomponent gaseous mixture, a numerical formulation is needed which 
applies to the equilibrium and near -equilibrium regions of the flow as well as 
the nonequilihrium region. The needed formulation is obtainable when a proper 
choice has been made for the dependent variables, and the rate equations in 
these variables are linearized over small intervals of time. By a properly 
chosen transformation, these rate equations can also be uncoupled, and the 
problem of obtaining expressions for the reaction vaxiables reduces to solving 
a system of linear, imcoupled, first-order differential equations. The solu- 
tion is, of course, easily obtained. However, in order to determine the 
needed transformation, it is necessary to solve an eigenvalue problem. It has 
been found that the needed transformation exists for all regions of the chan- 
nel, and a method has been developed which reduces the eigenvalue problem to 
finding the eigenvalues of a symmetric matrix. This greatly simplifies the 
computational problem because munerical methods for finding eigenvalues of a 
general matrix are unreliable in their accuracy. 

The integration procedure used in reference 10 can be regarded as a 
special application of the present method. In that work matrix methods were 
not required because of the low rank of the system. The method is presently 
being applied to a more general case. It appears that the necessary matrix 
manipulations can be performed accurately, and without adding greatly to the 
computing time per case experienced in reference 10. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., Oct. 28, 1965 
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APPENDIX A 


PROPERTIES OF THE MATRIX A 


The r X r matrix A = [aiz] of equation (36) has its elements defined 
in equations (7 ) • 


9-iZ - “ ^i)o 




7j / 


7-1 s ^Zj 

.n 7j 




J=i ^O 


(Al) 


This can he written 


^iZ “ Pi^iZ^Z 


where 


Pi = (1 - ^i)o^ 1 z = 


and 


1^iZ - 


Y PjJ^Zj 
L 7 j 




It is seen that h^^ = h^i and hence the matrix B = [hi^] is an r X r real 
symmetric matrix. Since Xi < it follows that > 0 (i - 1 ^ 2 ^ . . . ^ r) . 

Likewise^ from its definition^ > 0 (z 1^ 2 ^ . . - ^ r) . 

Let P - diag (pi^, P2, . . . ^ Pr) aJ^d Q diag (qi^ ^ q^) . Then^ 

from the construction of and it follows that 

A - PBQ 


Both P and Q are diagonal positive definite real matrices and B is real 
symmetric . 

Theorem 1 . 

Let A = PBQ where B is real symmetric and where P and Q are real^ 
diagonal^ and positive definite. 

P = diag (pi, P2^ • • Pr) ^ 

Q = diag (q.1, q^^ • • • ^ < 3 r) > 


Pj_ . «^r 

q.1 ^0^ i=-1^2^...,r 
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Then^ the eigenvalues of A are the same as the eigenvalues of a real symmet- 
ric matrix 

F DBD 

where 

D =■ diag Cj PiQ. 1 ^ s/PgQs ^ ^ %/Pr^x) 


Proof : 

Since P and Q are diagonal with rank r, each has an inverse. It is 
possible^ then^ to define a matrix 


Also^ let 
Then 

Now 

where 


s“^ = PQ“^ = diag (pi/qi, Pa/q^, • • • Pr/lr^ 

T = QPQ 

A = PBQ = PQ"^QBQ = S”^T 
lAl - A] = lAI - S-^T] = IS-i]tAS - T] 

S = diag (qi/pi^ Oa/Pa^ ■ • 


Consequently 


lAI - a| = |S“^| . 1?\S - t] 

From the definition of It Is seen that jS"^] ^ 0. Thus |?\I - A| = 0 If 

and only if |?\S - t| = 0. This shows that the eigenvalues of A are pre- 
cisely the roots of |?\S - T| =0. 

Now define the matrix R by 

R = diag C/pi/ai^ •Jvziqz, ■ ■ ■> s/Pn/<3n) 

This definition is made so that RSR - I. Consequently 

Ir] • |as - t[ . 1r| = |ai - rtr| 

Again, |r| 0 so that the roots of |AS - t| = 0 are the eigenvalues of the 

matrix 

F = RTR 

This shows that the eigenvalues of F are the same as the eigenvalues of A. 
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It -will be shown that F is symmetric. Since and B are all 

symmetric^ 

R^ = R^ = Q, and = B 

Then 

F^ = = R^Q^^Q%^ = RQBQR = F 

Since R and Q are diagonal^ they commute. Let 


D = RQ =- OR = C/piii, s/Ps'ls^ • • • ^ \/Pn1n) 


Then 


F = DBD 


Q.E.D. 


Theorem 2 . 

Let A and F be the matrices of theorem 1, and let R_^= diag (\/pi/q.i^ 

si Pg/qg^ • • • ^ s/p^/q.j,) . If A is an eigenvalue of F and t its correspond- 

ing eigenvector, then t] = R| is the eigenvector of A corresponding to A. 
Proof: 

— ► 

Theorem 1 shows that A is an eigenvalue of "both F and A. Since ^ 
is the corresponding eigenvector of F 

F| = (RQBRQ)T = A? 


Thus 


PQ"^R“^(RQBRQ)T = APQ"^R“^| 


The matrices R and Q commute since they are both diagonal. 


Hence 


(PBQ)r| = APQ"^R"^| 

Now PQ"iR"i = R, since PQ"^ = diag (pii/q.i, Vzt^> 
A = PBQ, so that the above equation is 


Define rf by 


AR| = AR? 

n - rT 


Pr/^r) • Also 


Then 


At) = At) 


so that if is the eigenvector of A corresponding to A. 


Q.E.D. 
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Theorem 3* 


The eigenvectors (i = 2, . . . ^ r) of the matrix A = PBQ of 

theorem 1 are linearly independent. 

Proof: Let F - DBD he the matrix used in theorem 1. Since F is real 

symmetric^ there is a nonsingular matrix H such that 

H'^ = dlag (Ai, A 2 , . . . , A^) 

where Ai, Ae^ • • • ^ Aj. are the elgen-valoies of the matrix F- The possihll- 
Ity that some eigenvalues may have multiplicity greater than unity is not 
excluded. Let L = dlag (Ax> Aa^ • • Aj.). Then 


FH = HL 

— 5 ^ — > — > 

Denote the column vectors of H by These column vectors 

of H are the eigenvectors of F^ for 

= Aili (l = 1, 2, . . r) 

— ^ 

It is to be noted that the are linearly independent since H is non- 

singular . 

Nov consider the matrix 

G - PH 


where 


dlag (v/pi/q. 1 , v/ps/cLa^ • • • ^ n/p^T^) 


The matrix G is non singular since it is the product of tvo nonsingular 
matrices, R and H. Thus, its column vectors form a linearly independent set. 
But the column vectors of G are 

T]i ~ 

The result of theorem 2 shows that these are the eigenvectors of A ^ PBQ. 

Q.E.D. 
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